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O ■ ABSTRACT 

(n: 

Vh ' We present a general linear algorithm for measuring the surface mass density 1 — k 

^ ' from the observable reduced shear g = 7/(1 — k) in the strong lensing regime. We 

show that in general, the observed polarization field can be decomposed into "electric" 
and "magnetic" components, which have independent and redundant solutions, but 
orthogonal noise properties. By combining these solutions, one can increase the signal- 
^ \ to-noise ratio by \/2- The solutions allow dynamic optimization of signal and noise, both 

^S) \ in real and Fourier space (using arbitrary smoothing windows). Boundary conditions 

^-^ ' have no effect on the reconstructions, apart from its effect on the signal-to-noise. Many 

f^ ■ existing reconstruction techniques are recovered as special cases of this framework. The 

11/ ■ magnetic solution has the added benefit of yielding the global and local parity of the 

Q^^ . reconstruction in a single step. 

^! Subject headings: (cosmology:) gravitational lensing 

^: 

c/3 . Gravitational lensing is rapidly providing large, independent data sets which are direct mea- 

• • , sures of gravitational potentials. One would like to understand the properties of gravitational 

.,-H i lensing reconstruction to optimize the solutions in a similar fashion as has been done for galaxy 

rS \ surveys (Vogeley and Szalay 1996). In the weak lensing regime, the reconstruction is linear, and 

d ' optimization of signal-to-noise becomes a straightforward problem. When lensing becomes strong, 

the equations appear to become non-linear, and the solutions are more difficult to understand. 
Of particular interest is the understanding of error propagation, to quantify the confidence of the 
solution, and to optimize the signal-to-noise for real, noisy data sets. The measurement of galaxy 
polarization, or reduced shear, yields two observable quantities. The ellipticity e = (1 — -R)/(l + R) 
is defined in terms of the minor/major axis ratio R, and the orientation of galaxies 93 is measured 
relative to the x-axis at each pixel of the reconstruction. One constructs a spin-2 vector with com- 
ponents ei = ecos2(/j, 62 = e sin 2^9. The reduced shear gi = Cj in the even parity case (i.e. for 
weak lensing), but changes to gi = Cj/e^ across a critical line. In the single source-lens plane lensing 
problem, the only unknown is the dimensionless surface mass density k. One would thus expect, 
in general, two independent solutions to exist, each using only half the observable information. 
If the noise properties in each solution are understood, and orthogonal, the two solutions can be 
combined in an optimal fashion. It is the purpose of this paper to systematically construct such 



solutions. We also quantify the existing solution algorithms in the same framework, and discuss 
error properties. 

We assume single source redshifts, which is a reasonable approximation for low redshift lenses. 
We will denote vectors and matrices by either bold symbols or explicitly using indices. One proceeds 
as follows: Let 

r.f^^ ^^ ] (1) 

y 72 -71 J 

and the observable reduced shear is given as G = r/(l — k), i.e. (1 — K)ga = 7a- There are two 
observables gi,g2, but only one independent unknown k. 

Let us first solve for 1 — k by noting that T has no magnetic component (Seljak 1997), i.e. 

Ti^ = {2did^V-^ - 6i^)d^djV-^Tij (2) 

which we can recast in terms of the observable reduced shear 

(1 - k)Gi^ = Ei^ = i2dld^V'^ - 5lm)^^^jV'\l - K)G^J. (3) 

We note that this is a linear differential equation for 1 — k in terms of the reduced shear. It is 
formally a fourth order equation. A peculiarity of equation (3) is that it is both overdetermined, 
and singular. In other words, 1 — k is the null eigenvector of the magnetic equation 

(1 - k)Bi^ = (1 - k)Gi^ - Eira = 0. (4) 

In the presence of noise, equation (3) may not have any solutions. We can choose instead a linear 
least squares problem where we minimize the noise S^: 

S^ = f BijB'^d^^. (5) 

S is a quadratic function in 1 — k, which we wish to minimize. Because of the global invariance 
transformation (Falco et al. 1985, Seitz and Schneider 1995), we cannot determine 1 — k up to 
a multiplicative constant, and in fact S would naively be minimized by k = 1. We impose a 
Lagrange multiplier L= J(l — K)^(i^x —1 = 0, and minimize S^ + AL for 1 — k and A respectively. 
The first variation leads to a linear eigenvector problem of the form A(l — «;) = A(l — k), while the 
The linear matrix A can be written as 






Kb ^ ....„ x.o.,.„ ^ (6) 



I Y, Giji^a? j 5(x, - x^) - 2G,,(x,)Gi„(x^) j e'^<^--^f^^hkjkik^d^h (7) 



where k^ = ki/\k\. In real space, the integral becomes 

—2 r 

+ -G(Xa)'^6{Xa - X/j) 



e^(xQ,)(5(xQ, - x^) + ^<^ [ei(xa)ei(x^) + e2(xQ,)e2(x^)] 

-2[ei(x«)cos(2^) -e2(x«)sin(20)][ei(x/3)cos(2^) -e2(x/3)sin(2^)] \ (8) 



where x = (xq, - x/3)/|xq, - x^| = {cos(0),sin(6')} and G = JY.ij Gfy 

A is formally infinite on the diagonal, but we can rescale it onto a new matrix which is zero 
on the diagonal: 

"^-^ = G(x„)G(x,) 2 • ('^ 

We note that the smallest eigenvector w of M is given by f = (1 — k,)G. One should note that across 
a critical line, the observable G changes parity, and one observes G/G^ instead (see for example 
Kaiser 1995, herafter K95). One of the main features of this method is its continuity across critical 
lines. We simply write down both solutions, 

Ki = 1 - v/G 

K2 = l-v/G^. (10) 

In a reasonably large field of view, where one knows the edges to be outside the outer critical line, 
one uses ki. At the point where ki = K2, we know that we have encountered a critical curve, and 
we switch the solution variable to K2- We will need to do the same across the inner critical curve. 
The local parity ambiguity can now be solved after the global solution, which is one of its main 
attractions. We will discuss a local deterministic procedure to determine the parity of (10) below. 

A direct solution for the smallest eigenvector of M can be computationally expensive. For 
an image which is N pixels on each side, 0{N^) operations are required for direct matrix solvers. 
This is the same operation count as directly solving the Kaiser and Squires (1993) (hereafter KS) 
procedure without the use of a fast Fourier transform (FFT), which becomes necessary if non- 
periodic boundary conditions are imposed. Linear algebra is, fortunately, a well exploited subject, 
and highly optimized and parallel libraries are available to find the required eigenvectors for matrices 
as big as 2^^ on a side in a day using a modern 100 CPU multiprocessor, which is sufficient to directly 
reconstruct images with 256^ pixels. In any case, it would clearly be desirable to use an iterative 
method, where each iteration would only involve a convolution. Since we are interested in the 
smallest eigenvector, an inverse power method would yield rapid convergence. Each such iteration 
involves solving a linear system, which is again straightforward using an iteration, since we know 
A to be positive semi-definite. The actual iterations are in fact just convolutions, which could also 



be accelerated using FFTs. The convergence rate for the inverse power method is given by the 
ratio of second smallest eigenvalue to smallest eigenvalue A2/A1, which approaches infinite speed 
for small noise. 

The standard "electric" mode reconstruction can be applied in a similar fashion, but requires 
prior knowledge of the local parity. Let u = 1 — k, Hij = Gij + 5ij, then the lensing equation reads 
(in the even parity case) 

didjuHij = (11) 

for which we can again define a least squares action S = f (\7^'^ didjuHij)'^ , resulting in a matrix 



G2(Xn 



^a/3 = — {-2x'x^ [Giji^a) + Gy(x^)] + Gy(x„)G/m(x/3) 5udjm - Ax'x^Sjm + 4x*x^x'x 

+(5(Xo - X/3 



1 + 



(12) 



where r = |xq, — x^| . We note the following property of minimum eigenvector solutions to (12): if u 
is a solution for a given Hij, then it is also a solution for addition of magnetic noise, Hij + N^- /u. We 
are considering noise fields N which are arbitrary traceless tensor fields, which can be decomposed 
into electric and magnetic components, in analogy with Equations (3,4). This is the opposite 
property of the magnetic solution for (9), where one could add electric noise to Gij + Nf- /v and 
keep any solution v invariant. We see that the two solutions have orthogonal dependencies on the 
noise, and thus expect their combination to improve signal-to-noise ratios by \/2- If the noise is 
known to contribute equally to E and B, as most sources would be expected to, we can use the 
difference between the two solutions as an estimate of the noise. 

An elegant observation by Kaiser (K95) was the realization that djuHij is a curl-free vector 
for any true solution n, allowing an integration of (11): 

dk\og{u) = -H-ld,Hij. (13) 

We can also solve an equivalent least squares problem by setting S = J \7^'^ {djuHij)[di^uHij^)(P'K. 
This strategy was explored by Lombardi and Bertin (1999), who considered accelerated direct 
solutions. The matrix is 

^a/3 = ^{x'x^[Gij{Xa) + Gij{Xf3)+Gii{Xa)Gij{xp)]} 



+6{Xa - X^ 



1 + 



G(xq, 



\2 



(14) 



One also obtains A''^ by a contraction of the appropriate term in A-^ (12). As K95 pointed out, 
the parity can be directly determined from the curl of a vector. Let us define I'i = djuGij. We can, 
in analogy to the magnetic solution, solve for both the parity and k by requiring that V x 1/ = 0. 
The corresponding action is S**^ = J(V^^V x v)"^, which results in a linear combination of the 
electric and Kaiser matrices, A*^ = A^ — A^. K95 proposed a similar procedure to determine the 
parity, which is necessary to construct equation (13) from the observable ellipticities Cj. 



It is instructive to examine the weak lensing limit. When G ^ 1, the eigenvector of Hij 

can be considered as a perturbation on the background H^- = 5ij with H^- = Gij. We note that 

Equations (12) and (14) agree to 0{G). We will choose a Fourier weight function k'^ in the action, 

so S^ = f{diu)d^u, and S^ = J V~'^{didjuGij)'^ . This breaks the degeneracy of eigenvectors in A'^, 

giving us A^ = —5"{r). The eigenvectors are the Fourier modes exp(ik-x) with eigenvalue k"^. The 

zero eigenvector is given by u^ = 1. The perturbed lowest eigenvalue is f kikjGij (see for example 

Schiff 1971, section 31). The perturbed zero eigenvector is given in Fourier space as the matrix 

element 

r A^ 

u^{ka) = / exp(-ikQ, • :s.a)-^d^^i3d^^a (15) 

We then have u^ = kikjGij/k"^, which is exactly the KS solution. The B equation (8) does not have 
a weak lensing limit, since all entries are 0{G'^). 

One can combine the algorithms by adding the actions, for example using S = S^ + S^ or 
any linear combination of the three actions described above. The resulting matrix A just becomes 
the sum of the matrices. An interesting combination is S = S + 25 , which has no quadratic 
terms in G except on the diagonal. We note in passing that this S'^ = "^Ziji^S)"^ i^ '^^ ^^^^ ^^^ least 
squares action, or likelihood, of the appropriately weighted reconstructed reduced shear as used in 
"maximum likelihood reconstruction" (Bartelmann et al. 1996): 

Sf^ = V-2 Ud, - ^-fvA K-{1- k)Q, = 0. (16) 

At each point, S'~^ /{I — k)^ is just the difference squared between the observed reduced shear Gij 
and the reconstruction for a given k field. 

The construction of S requires no prior knowledge of parity, while S does, so one might 
expect to proceed by first solving S^ and then using the weak lensing equation directly. The B-type 
solution (9) has only made use of the B component of the reduced shear, which works well in the 
strong lensing regime, but results in very poor optical depth estimation in the weak lensing regime. 
One can apply a linear method, such as KS on the reconstructed (unreduced) shear 

r° = f. (IT) 

No parity knowledge is needed for equation (17), which automatically returns the correct parity for 
r. We can then reconstruct k^^ from F, and compare this second solution with (10). This allows 
us to locally and globally determine the parity of the solution. In the final map, one can combine 
the KS inferred k with the k° inferred using the new procedure. 

The construction of M from A can be generalized to optimize for prior knowledge of noise and 
signal statistics. Equation (5) can be generalized with arbitrary filter weightings in Fourier space, 

S= f Bij{k)B'^-'k)W{\k\)d^k. (18) 



The window function W{k) transforms to W{r) in real space. The matrix entries are accordingly 
modified as follows: 

A^^ = G,,(x„)Q,(x^)[H^(r) + 2yi(r)-2V2(r)] 

-2Gij{^^)GU^p)Ux'x'6,„. [6Fi(r) - 3V2{r) + V^ir)] 

+x'x^x^x"' [-15Vi(r) + 15V2{r) - 6V3(r) + Viir)] 

+G,,(x,)G,„(x^)|4£*x'(5,„ [6Vi(r) - 3V2(r) + V^(t)\ 
+x^x^x'x™ [-15Vi(r) + \hV2(r) - ^Vz{r) + V^{r)\ 

X {5ijW2{r) + x'x^ [W{r) - 2W2{r)]} (19) 

where Vn{r) = r^~^d^\7~^W{r) and Wn{r) = r~"' J^ s^~^W{s)ds. The original formulae (8,12,14) 
are reproduced for W{r) = 6{r) (2-dimensional). In addition, we can impose a weighting in real 
space, as done in Equation (9). Instead of dividing to obtain a unit diagonal, one might also 
envision weighting by the expected local signal-to-noise (Pen 1999). In general, we see that setting 
Maf3 = Aaf3'w{xa)w{'x.i3) relates the zero eigenvector u(xq) of Aafj to the computed zero eigenvector 
m{xa) of Maf3 by u(xa) = m{xa)'w(xa). We see immediately that no boundary condition artifacts 
are ever introduced when one truncates A at arbitrary boundaries, for example by choosing w to 
be 1 on the domain with data and elsewhere. One should note that if too much noise is added 
to S, it can happen that the smallest eigenvector is dominated by noise, while the second smallest 
eigenvector actually contains the correct solution. 

We have presented a direct linear solution strategy for the strong lensing problem. The B- 
type solutions work contiguously across critical lines and critical density lines without prior parity 
knowledge. The critical lines are then obtained from the solutions (10) themselves. The existence 
and uniqueness of the solutions and the critical lines are therefore clearly established, which is 
in general not clear for non-linear solution schemes. Its noise properties depend only on the B- 
type noise, in contrast to the KS weak lensing and the E-type reconstruction procedure which 
depends only on E-type noise. We have fully exploited the fact that the two observables gi,g2 allow 
the construction of two independent solutions with orthogonal noise properties. Combining them 
explicitly allows one to improve the signal-to-noise by up to \/2. In addition, it allows for arbitrary 
weightings in both real and Fourier space. We have shown how to apply the same procedure to 
the K95 and standard maximum likelihood algorithms. The solutions furthermore allow arbitrary 
weightings in both real and Fourier space, which allows optimization of signal-to-noise in the final 
reconstruction. 
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